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ABSTRACT 


The  Terrain  Shielding  Algorithm 
Used  in  The  MASR  Airborne  Surveillance  Radar  Program 


The  algorithm  described  herein  quickly  determines  if  a  point  on  the 
surface  of  the  earth  is  visible  or  not  from  a  remote  aircraft  due  to  masking 
by  nearby  terrain  features.  It  was  designed  to  work  in  an  airborne  radar 
surveillance  system,  updating  a  16  x  16  km  area  every  4  seconds  with  a 
resolution  of  2S0  x  250  m. 
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Z.  INTRODUCTION 

The  MASR1  was  am  agile  beam  standoff  airborne  surveillance  radar  (ASR) 

system  that  used  a  wide-band  data  link  to  send  target  reports  to  a  Ground 
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Processing  Center  (GPC)  for  tracking  and  display.  Typically,  the  MASR  would 
be  flying  at  an  altitude  of  10,000  feet  30  km  from  the  target  and  would  be 
examining  several  independent  2  x  2  km  Areas  of  interest  (AOI)  that  could 
be  placed  anywhere  inside  a  fixed  16  x  16  km  area  by  using  an  interactive 
graphics  display 

In  such  a  situation,  it  is  desirable  to  know  whether  the  AOI's  are 
visible  from  the  aircraft  so  that  one  would  not  waste  the  ASR's  resources 
examining  shielded  areas,  and  to  optimize  the  flight  path  in  the  mission  planning 
stage,  or  to  anticipate  track  dropouts  and/or  how  long  a  track  will  coast.  An 
AOI  can  become  invisible  because  it  is  in  a  valley,  or  2m  intervening  mountain 
masks  it,  or  it  is  shielded  by  trees,  foliage,  buildings,  etc. 

In  order  to  calculate  the  terrain  shielding,  one  must  have  a  hipsographic 
data  base.  Fortunately,  we  were  able  to  obtain  from  the  Defense  Mapping 
Agency's  Topographic  Center  in  Washington,  DC,  a  magnetic  tape  containing 
the  elevation  of  Eastern  Massachusetts  for  an  array  of  points  spaced  63.5 
meters  apart.  These  data  include  only  the  hipsographic  (elevation)  data, 
nothing  for  trees  or  cultural  details,  but  these  could  be  added  to  the  data 
base  if  desired.  Appendix  1  describes  in  more  detail  the  character  of  this 
hipsographic  data  base. 

Since  the  aircraft  is  moving,  the  terrain  that  is  shielded  changes  as  the 
aspect  from  the  aircraft  changes,  consequently,  a  real-time  computation  of 
this  effect  is  required.  With  the  MASR  aircraft  flying  at  speeds  of  about 
100  knots  over  Central  Massachusetts,  the  terrain  that  is  shielded  remains 
essentially  constant  for  a  few  tens  of  seconds  using  an  elemental  area  (or 
precision)  of  250  x  250  m.  This  time  constant  would  be  less  for  more 
mountainous  terrain,  greater  for  flatter  terrain. 


One  of  the  functions  of  the  GPC  was  to  provide  this  terrain  shielding 

information  in  real  time.  This  was  done  by  devoting  a  background  task  to 

this  job.  The  GPC,  using  a  Fortran  task  in  a  DEC  PDP  11/70  computer,  was 
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able  to  refresh  a  16  x  16  km  display  in  4  seconds  with  250  m  resolution. 

Since  this  was  fast  enough  for  the  MASK  requirements,  no  effort  was  spent 
trying  to  speed  up  the  program  by  writing  it  in  machine  language  etc. 

The  resolution  and/or  the  total  area  could  have  been  increased  at  the 
expense  of  more  storage  for  the  hipsographic  data.  The  PDP  11,  being  a  16  bit 
computer  with  32  k  words  of  address  space,  is  hard  pressed  to  increase  the 
core  resident  hipsographic  data  beyond  the  16  k  words  that  were  used.  A 
larger  hipsographic  data  base  could  be  stored  on  a  disk  and  selective  portions 
read  in  as  required ,  but  this  would  complicate  the  program  and  slow  it  down 
so  it  would  not  be  suitable  for  our  real  time  application.  A  32  bit  computer 
would  ameliorate  this  problem. 

11.  ALGORITHMS 

A.  Basic  Algorithms 
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The  basic  algorithm  is  a  ray  tracking  approach  developed  by  RADC  .  It 

starts  with  a  point  on  the  surface  of  the  earth  and  constructs  a  straight 

line  (ray)  to  the  aircraft.  The  program  "moves  up"  this  ray  in  constant 

increments  (200  meters)**.  At  each  increment,  the  ray  altitude  is  computed  and 

compared  with  altitude  of  the  surface  directly  beneath  it.  If  the  surface 

altitude  exceeds  the  ray's  altitude,  the  initial  point  is  considered  blocked. 

If  not,  the  next  point  up  the  ray  is  examined, etc. ,  until  the  ray  becomes 

higher  than  the  highest  point  in  the  hipsographic  data  base.  In  this  case, 

the  initial  point  is  considered  to  be  in  the  clear.  The  program  then  moves 

*• 

along  the  earth's  surface  a  constant  increment  (250  meters)  to  start  a  new 
initial  point  and  the  process  is  repeated.  When  the  entire  AOI  has  been 
scanned,  250  meters  at  a  time,  the  algorithm  terminates. 

•The  time  required  is  directly  proportional  to  the  area  and  inversely  pro¬ 
portional  to  the  square  of  the  resolution. 

••These  constants  were  chosen  to  be  commensurate  with  the  resolution  of  the 
hipsographic  data  base. 
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B.  RADC  Implementation 

The  code  that  RADC  developed  to  inclement  this  algorithm  started  with 
geodetic  spherical  coordinates  (longitude,  latitude,  altitude  above  mean 
sea  level  -  X,  8,  h)  for  the  aircraft  and  first  surface  point.  These  would 
be  converted  to  a  geocentric  rectangular  coordinate  system  (x,  y,  z  from  the 
center  of  the  earth)  so  that  the  ray  and  increments  along  the  ray  are  trivial 
to  calculate.  For  each  increment  up  the  ray,  the  x,  y,  z  position  would  be 
converted  to  X,  8,  h  and  the  resultant  "h"  would  be  compared  first  to  the 
maximum,  then  to  the  surface  altitude.  This  second  coordinate  conversion  was 
necessary  because  the  hipsographic  data  were  stored  as  of  function  of  X,  8. 

After  the  ray  is  shown  to  be  blocked  or  in  the  clear,  the  next  surface 
point  (X,  8)  is  found  (by  adding  a  constant  AX  and  A8  converted  to  x,  y,  z  and 
the  process  repeated  until  the  desired  areas  has  been  scanned.  The  two  con¬ 
versions  from  X,  8  to  x  y  z  and  back  again  X,  8,  were  done  to  take  the  curva¬ 
ture  to  the  earth  into  account. 

4096  surface  points  are  required  to  scan  16  x  16  km  with  a  precision  of 
250  meters;  each  surface  point  required  one  conversion  from  spherical  to 
rectangular  coordinates  and  generated  one  ray.  On  the  average,  there  might 

4 

be  10  increments  along  each  ray  to  be  examined,  so  some  4  x  10  conversions 
from  rectangular  to  spherical  coordinates  were  required  in  all. 

The  conversion  from  rectangular  to  spherical  coordinates  (requiring  at 
least  2  arc  tangents  and  2  square  roots)  is  clearly  the  limiting  factor  in 
this  implementation  for  any  real  time  application. 

C.  MASR  implementation 

The  MASR  Implementation  of  the  terrain  shielding  algorithm  substitutes 
for  the  coordinate  conversion  an  expression  that  is  much  faster  to  calculate. 

In  addition,  it  not  only  stays  in  one  coordinate  system  for  all  the  terrain 
shielding  calculations,  but  keeps  it  throughout  all  aspects  of  the  MASR  soft¬ 
ware.  The  output  display  showing  roads,  political  boundaries,  etc.,  the  position 
of  the  aircraft  and  its  flight  path,  and  the  targets  and  beam  prints  were  all 
shown  on  the  same  projection. 
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The  obvious  choice  for  this  projection  is  a  gnomic  projection  to  a  plane 
tangent  to  the  earth  at  the  center  of  the  16  x  16  km  area.  The  chief  attrib¬ 
ute  of  this  projection  is  that  all  straight  lines  are  great  circles  and  vice 
versa,  so  that  one  can  follow  the  path  a  radio  wave  takes  by  adding  a  constant 

A  x,  A  y  to  the  starting  point. 

* 

The  gnomic  projection  obviates  the  need  for  back  and  forth  coordinate 
conversions  for  the  calculations  involving  positions  on  the  ray  in  the  x-y  plane 
However,  the  curvature  of  the  earth  introduces  a  complication  in  determining  the 
height  of  the  ray  above  mean  sea  level.  Figure  II.  1  shows  a  plane  passing 
through  the  aircraft,  center  of  the  earth,  and  the  starting  point  of  the  ray 
on  the  ground.  From  the  law  of  cosines: 

(a  +  r)2  »  d2  +  (r  +  h  )2  -  2d  (r  +  h  )  cos  o 

o  o 

(h  +  r)2  «  y2  +  (r  +  h  )2  -  2y  (r  +  h  )  cos  u 

o  o 

where-,  a  is  the  altitude  of  the  aircraft; 

hQ  is  the  altitude  of  the  terrain  at  the  rays  starting  point; 

h  is  the  altitude  of  the  ray  a  distance  of  y  from  the  starting 
point; 

d  is  the  distance  along  the  ray  from  the  aircraft  to  the  start¬ 
ing  point; 

a  is  the  angle  APO; 

r  is  the  local  radius  of  the  earth. 

*For  small  enou^i  areas,  where  the  difference  between  projection  is  less 
than  the  precision  of  the  hipsographic  data  one  can  use  a  projection  that 
is  more  convenient  from  other  respects.  In  MASR,  for  example ,  the  largest 
area  of  concern  for  terrain  shielding  was  16  x  16  km  and  the  hipsographic 
precision  (quantization)  was  125  meters.  Consequently,  we  used  a  modified 
cylindrical  equal-spaced  projection  (linear  in  latitude  and  longitude)  since 
the  output  of  the  program  that  converted  distance  and  angle  from  the  aircraft 
to  ground  coordinates  was  latitude  and  longitude. 4  If  the  output  were  in 
rectangular  cartesian  coordinates,  one  would  use  an  orthographic  projection, 
or  if  military  use  demanded,  one  could  use  universal  Transverse  Mercator. 
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Simplifying,  this  becomes: 

(h  +  r)2  =  y2  +  {r  +  h  )2  +  Ay 

o 

where : 

A  =  (a  +  r)2  -  d2  -  (r  +  ty>)2 
d 

Solving  for  h,  we  have: 
h  =  -r  +  \/  (r  +  hQ) 2  +  y2  +  Ay 

“  +  h° + 1/2  ^ 

Dropping  the  last  term  causes  error  less  than  1  foot,  so  we  have: 

h  =  h  +  x2  -t 

0  2(r  +  h  ) 

o 

As  y  goes  from  0  to  d,  h  goes  quadractically  from  h  to  a.  If  the  ray's 
starting  point  is  beyond  the  horizon  (a  <  90°) ,  A  <  0  and  h  will  start  off  de¬ 
creasing  as  we  move  up  along  the  ray.  Such  a  point  will  be  blocked  unless  the 
terrain  altitude  falls  off  at  a  greater  slope. 

An  approximate  correction  for  the  curvature  of  the  ray  due  to  refraction 
can  be  made  by  making  r  larger  than  the  actual  earth's  radius  by  a  factor  of 
1.0  to  1.5  depending  on  the  actual  index  of  refraction  (4/3  is  a  very  good 
nominal  value  to  use) .  More  accurate  corrections  for  refraction  or  the  ellip- 
ticity  of  the  earth  are  not  worthwhile  because  the  hipsographic  data  base  is  so 
coarsely  quantized. 

In  effect,  we  have  replaced  the  conversion  from  rectangular  to  spherical 
coordinates  by  the  much  simpler  quadratic  expression  given  above. 

D.  Other  Implementations 

The  other  two  methods  of  implementing  terrain  shielding  will  be 
briefly  discussed.  Neither  of  these  methods  was  used  in  the  MASR  experiments. 
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1.  Perspective  Projection 

This  algorithm^ constructs  a  continuous  three  dimensional  function 
from  the  hipsographic  data  which  is  viewed  from  the  aspect  of  the  aircraft. 

A  hidden  line  algorithm  marks  the  area  that  is  shielded.  In  addition  to 
giving  a  visible-invisible  plot,  this  method  also  gives  a  perspective 
projection  of  the  topographic  data,  which  might  be  very  useful  in  the 
mission  planning  stages.  Reference  7  gives  a  detailed  discussion  of  this 
method. 

2 .  Fixed  Scenes 

This  scheme  uses  one  of  the  above  methods  to  generate  shadow  plots  for 
several  candidate  aircraft  positions.  These  plots  are  discretely  switched 
as  the  aircraft  moves  about.  This  method  is  satisfactory  for  fixed  flight 
paths  and/or  not  too  hilly  terrain.  Potentially,  this  is  the  fastest  method 
for  real  time  use  from  a  computational  point  of  view  (presumably  a  library 
of  fixed  scenes  would  be  prepared  and  stored  on  a  disk  and  a  simple  look  up 
algorithm  based  on  the  aircraft’s  location  would  requisition  the  closest  one), 
but  it  lacks  generality  and  is  not  an  aesthetically  satisfactory  solution. 

III.  THE  TERRAIN  SHIELDING  PROGRAM 

A  tutorial  version  of  the  terrain  sheilding  program  has  been  prepared 
with  the  interactive  graphics  display  removed  so  that  only  the  essentials 
remain  .  A  flow  chart  of  this  program  is  shewn  in  Figure  III.l.  A  complete 
listing  is  given  in  Appendix  2.  All  that  is  required  to  make  it  a  working 
program  on  any  installation  is  to  have  a  subroutine  to  fill  up  the  MAP  array 
with  the  hipsographic  data  in  the  format  described  in  Appendix  1.  This  would 
go  in  place  of  statements  {9(129  to  0(131. 

The  innermost  loop  (statements  0071  to  0080)  is  the  area  where  the  pro¬ 
gram  spends  most  of  its  time,  so  the  code  has  been  arranged  to  try  to  speed 
up  the  calculations  in  this  area.  That  is  why  there  is  an  "odd"  system  of 
units  used  —  height  inoontour  intervals,  and  distance  in  units  of  the  scope 
display  instead  of  meters  or  feet. 
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Two  examples  of  the  output  of  this  program  are  shown  in  Figures  III. 2 
and  XII. 3.  Both  figures  show  the  terrain  that  is  shielded  (as  an  asterisk) 
from  an  aircraft  50  km  due  east  of  the  center  with  a  depression  angle  of 
2.3°.  The  resolution  is  250  m  in  Figure  III. 2  and  125  km  in  Figure  III. 3. 

The  center  of  the  figure  corresponds  to  42°  30'  N,  71°  37.5*  W.  West  is  up, 
north  is  the  right  due  to  the  way  the  scanning  raster  was  set  up  in  the 
program. 

In  Figure  III.  4  the  aircraft  has  been  moved  1  km  north  of  its  position 
in  Figure  III. 3.  This  position  is  where  the  MASR  aircraft  would  be  16 
seconds  later  if  flying  at  120  knots.  Sixteen  seconds  is  the  time  it  takes 
the  program  to  completely  calculate  the  shadowing  with  125  meters  resolution. 

A  close  examination  of  these  two  figures  reveals  essentially  no  difference 
showing  that  this  program  runs  fast  enough  for  real  time  application  to  the 
MASR  aircraft  in  the  Eastern  Massachusetts  environment.  This  satisfactory 
situation  is  made  even  better  (by  a  factor  of  4)  when  250  meters  resolution 
is  used. 

If  this  program  is  run  with  a  "ft"  as  the  logical  unit  for  the  output 
device,  there  will  be  no  shadow  map  generated.  The  program  was  used  in  this 
manner  to  generate  the  data  shown  in  Figure  III. 5  which  is  an  average  of 
the  percentage  visible  from  an  aircraft  flying  in  the  four  cardinal  directions . 
from  the  center. 
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APPENDIX  1 


Hipsographic  Data  Base 

The  Defense  Mapping  Agency  has  transformations  of  elevation  information 
from  maps  or  photos  to  x,  y,  z  digital  coordinates  available  on  1/2  inch 
magnetic  tape.  Approximately  two  2400  ft  reels  of  tape  are  required  to  record 
the  elevation  data  corresponding  to  one  1:250,000  scale  topographic  map  (ten 
thousand  elevation  points  are  produced  per  square  inch) . 

The  elevation  data  are  derived  from  topographic  maps  by  manually  tracing 
contour  lines  on  a  Digital  Graphics  Recorder  (DGR)  which  digitizes  points  at 
which  the  contour  line  crosses  superimposed  grid  lines.  The  grid  system  is 
an  integral  part  of  the  DGR  and  consists  of  lines  etched  on  a  glass  plate  at 
intervals  of  0.01  inch  in  both  x  and  y  axes.  The  result  of  the  DGR  operation 
is  a  series  of  coordinates  which  define  a  contour  line.  After  arranging  the 
data  into  an  x-y  sequence,  the  next  step  is  known  as  the  Planar  process,  a  name 
derived  from  the  basic  procedure  which  uses  points  on  contour  lines  to  define  a 
plane  on  which  the  elevation  of  a  point  between  contours  must  fall.  The  end 
result  is  an  elevation  at  each  0.01-inch  grid  point  from  edge  to  edge  of  the 
source  map. 

The  resolution  of  the  elevation  data  has  an  absolute  and  a  relative  aspect. 
The  resolution  of  the  device  which  converts  line  (graphics)  information  to 
digital  forms  is  0.01  absolute  inch.  The  resolution  of  the  "absolute"  inform¬ 
ation  is  relative  to  the  scale  of  the  source  map.  For  example,  a  digital 
record  of  the  elevations  obtained  from  a  1:250,000  scale  map  would  have  an 
absolute  resolution  of  0.01  inch  which  is  a  relative  resolution  of  about  63.5 
meters  on  the  ground  because  of  the  source  map  scale.  The  accuracy  of  the 
elevation  data  is  purely  relative  to  the  accuracy  of  the  source.  Data  obtained 
from  a  Class  A  1:250,000  scale  topographic  map  are  neither  more  nor  less 
accurate  than  the  data  which  could  be  obtained  by  ceureful  scaling  and  inter¬ 
pretation  of  the  reproduction  material  which  has  not  been  distorted  by  the 
printing  process.  Also,  the  computer  process  which  interpolates  the  normal 
contours  on  a  topographich  map  to  produce  a  matrix  of  equally  spaced  elevations 
is  more  consistent  than  a  manual  process. 
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The  elevation  data  are  organised  in  an  x-y  matrix  with  only  the  elevation 
value,  z  recorded.  The  x  and  y  coordinates  are  implied  by  the  position  of 
the  z  within  the  record.  The  relation  of  the  x-y  coordinates  to  geographic 
position  is  included  in  the  file  for  each  map  sheet.  The  horizontal  and 
vertical  datums  are  the  same  as  those  on  the  source  map. 

The  first  step  in  using  these  data  is  to  convert  the  data  from  x-y 
coordinates  to  latitude-longitude.  This  was  done  by  using  the  two  expressions 

9*a  +  bx  +  cy  +  dxy 

X*e  +  fx  +  gy  +  hxy 

These  expressions  account  for  the  skew  in  the  original  topographic  map 
(a  Transverse  Mercator  projection  for  the  Boston  map) . 

The  first  file  on  the  tape  contains  the  x-y  coordinates  and  corresponding 
latitude  and  longitude  for  the  four  comers  of  the  area  that  is  digitized  on 
on  the  tape  (a  1°  by  1*  region) .  These  can  be  substituted  in  the  above 
equations  to  give  8  equations  with  8  unknowns,  a  through  h.  These  can  easily 
be  solved.  For  the  western  half  of  the  Boston  map,  these  constants  were: 

a  -  41.98486136 
b  -  3.062634379  x  lO-6 
c  -  5.710621830  x  lO-4 
d  -  2,848018909  x  10-11 
e  -  -72.07252640 
f  -  7.660974315  x  10“4 
g  -  -1.002041361  x  10“5 
h  -  7.157438295  x  lO-9 

The  second  step  is  to  get  the  desired  data  off  the  tape  that  covers  the 
area  of  interest.  These  data  can  be  converted  from  x-y  to  X-6  by  the  above 
relations. 


The  thirl  step  is  to  convert  these  data  back  into  an  x'-y'  using  the 
projection  that  the  terrain  shielding  algorithm  is  to  use  -  e.g.,  Gnomic. 

Lastly,  these  data  are  scaled  to  be  the  desired  resolution  (in  the 
MASR  experiment,  every  other  point  was  decimated  to  give  125  meters  resolution- 
a  compromise  between  the  storage  requirements  and  execution  time  vs.  precision) 
and  put  in  a  two  dimensional  Fortran  array  MAP  (I,J)  where  MAP  (0,0)  is  the 
southwest  comer's  elevation,  MAP  (NMAP,0)  is  the  northwest  comer's  elevation 
(MAP(NMAP,  nmap)  is  the  dimension  of  the  array  MAP) . 

In  the  MASR  program,  this  array  was  written  onto  the  disk  and  a  special 
subroutine  (using  QIO)  was  used  to  retrieve  it  directly  and  quickly  without 
going  through  the  DEC  disk  files  management  software. 

Figure  A. 1  shows  the  output  of  the  Hipsographic  data  base  for  20  x  20  km 
centered  on  42°30'n,  71°41.25'W.  The  density  becomes  darker  as  the  terrain 
altitude  exceeds  three  fixed  thresholds.  This  figure  illustrates  the  type  of 
terrain  used  to  generate  the  data  for  this  report.  The  altitude  varied  from 
200  to  600  ft  for  the  MASR  AQI. 

Although  the  data  base  contains  no  foliage  information,  provision  was 

made  in  the  program  to  include  a  very  rough  correction  for  the  effect  of 

* 

foliage.  The  approximation  was  that  all  trees  are  a  constitn t  height  ,  and 
uniformly  surround  every  ray- starting  point.  In  other  words,  the  starting 
altitude  is  decreased  by  the  height  of  the  trees  (line  0059  in  Appendix  2) . 

For  the  MASR  AOI,  a  densely  forested  area,  this  approximation  was  fairly 
good  except  for  points  on  an  interstate  highway  when  the  aircraft  was  di¬ 
rectly  in  line  with  the  road.  Unfomately,  this  was  the  situation  that 
existed  in  the  MASR  experiments,  so  this  feature  was  not  used. 


* 

An  input  parameter  (line  0024  in  Appendix  2) .  The  value  "0"  was  used  for 
the  MASR  experiments. 
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T.FTH  /UR 
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PROGRAM  TS  1 


LIVES  IN  UFD:  146,341 


TERRAIN  SHIELDING  PROGRAM 
BV  P.  L.  FLECK 


PARAMETER  NMAP-128 

DIMENSION  MAP<  NMAP ,NMAP ) ,MARK( 130) , IMARK(  65 ) 
I NTEGER*2  I MAPI NMAP+NMAP ) 

1 NTEGERM  BLKL4.MPBVTS 
BYTE  I FORK, MARK 

EQUIVALENCE  ( MAPI  1 , 1 > , IMAP  I  1 ) ) 

EQUIVALENCE  I MARKI 1 ) , I MARK! 1 ) > 


DATA  RE/8491 253 . / 
DATA  SIZMAP/16000./ 


{RADIUS  OF  EARTH  AT  42.5  DEGREES  •  4/3 
I  METERS  WIDE 

{THERE  ARE  2  CONTIGUOUS  .OLB  FILES 


/**»**  DEFINITIONS  FOR  SOME  CONSTANTS 


ALPHA  - 
BETA  — 
10,0) 
1X1 ,IY1 
MUST 
NMAP  -- 
RE  — 

SIZMAP 
I NMAP/ 2 
10,0)  I 


-  LATITUDE  OF  TANGENT  PLANE  FOR  GNOMOHIC  PROJECTION 
LONGITUDE  OF  TANGENT  PLANE  FOR  GNOMONIC  PROJECTION 
IN  SCOPE  UNITS  IS  THE  POINT  I  ALPHA, BETA) 

—  LOWER  LEFT  CORNER  OF  SCOPE  DISPLAY  IN  SCOPE  UNITS 
LIE  IN  THE  REGION  +  OR  -  4096 
DIMENSIONS  OF  MAP  ARRAY! =128 > 

RADIUS  OF  EARTH  AT  LATITUDE  ALPHA.  THIS  CAN  BE  MULTIPLIED  BY  4/3 
TO  APPROXIMATE  REFRACTION  EFFECTS 

—  EXTENT  OF  MAP  ARRAY  IN  METERS 

, NMAP/2  >  IS  THE  POINT  {ALPHA, BETA) 

S  THE  SOUTH-WEST  CORNER,  SIZMAP-0.707  FROM  {ALPHA, BETA) 


DEFINITIONS  FOR  CALCULATED  VARIABLES 


CONTUR  --  NUMBER  OF  METERS  PER  CONTUR  UNIT  I  PROBABLY  I  FOOT)  (-0.3048) 

TREE  —  HEIGHT  OF  FOLIAGE  ABOVE  TERRAIN,  IN  CONTOUR  UNITS 
XUTOM  —  CONVERSION  FACTOR  FROM  SCOPE  UNITS  TO  METERS 
UMAP  --  ONE  MAP  UNIT  IN  METERSI-125) 

XMAPS  —  THE  RATIO  OF  SCOPE  UNITS  TO  MAP  UNITS. I -1/16 ) 

RE  —  RADIUS  OF  EARTH  AT  THE  GIVEN  LATITUDE  ALPHA  I  IN  METERS) 

RES  --  RESOLUTION  I  IN  METERS)  OF  TERRAIN  SHIELDING  ON  THE  GROUND<-250  OR  125) 
XINTVL  --  RESOLUTION  ALONG  GROUND  IN  SCOPE  UNITS  132  OR  16) 

RAYINC  --  AMOUNT  WE  MOVE  ALONG  THE  RAY  FROM  THE  GROUND  TO  THE  PLANE  (-200) 

(IN  METERS) 

RINSCU  —  RAYINC  IN  TERMS  OF  SCOPE  UNITS  (25.6) 

SCALE  —  SCALING  FACTOR  FROM  LAT-IRAD)  TO  SCOPE  UNITS  -RE/XUTOM 
SUTOM  —  CONVERSION  FROM  SCOPE  UNITS  TO  METERS  (8000/1024-7.8125) 

SUTOM-8000. 0/1024.0 

TYPE  INPUT  LOGICAL  UNIT  NO  OF  OUTPUT  DEVICE  FOR  SHADOW  MAP* 

ACCEPT  *,LUN 
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FORTRAN  1V-PLUS  V02-51 
T.FTN  /WR 

C  IF  LUN-0,  NO  “SHADOW  MAP'  IS  PRODUCED 

C 

C  DETERMINE  RESOLUTION  TO  BE  USED 

CANT  BE  ANY  FINER  THAN  THE  INPUT  HIPSOGRAPHIC  DATA  BASE 
0014  1  TYPE  AND  RESOLUTION  TO  BE  USED  (1-125  OR  2-250  METERS)' 

0015  ACCEPT  *  , ORES 

0016  1 F ( 3RES.LE .0 .OR. ORES .GT.2 )  GO  TO  1 

0017  RES“125.0*ORES 

C 

0018  1X1  —  1024 

0019  IV 1  —  1024  IDO  ENTIRE  16X15  KM  ARRAY  TOO 

C 

0020  XINTVL-RES/SUTOM  1-32  OR  16 

0021  IDEL-XINTVL  ISCAN  INCREMENT  IN  SCOPE  UNITS 

0022  NINV-128/JRES  1NO  OF  POINTS  TO  EXAMINE  EACH  SCAN  (64  OR  128) 

C  AND  ALSO  THE  NO  OF  SCANS  AS  WE  ARE  DOING  A  SQUARE  ARRAY 

C 
C 

0023  CONTUR  -  0.3048 

0024  TREE  -  0.  /  CONTUR 

C  USE  0  METERS  FOR  TREE  HEIGHT 
0025  UMAP-SIZMAP/NMAP  1-125 

0026  XMAPS  =»  SUTOM  /  UMAP  1-1/16 

0027  RAYINC-200. 

0028  RINCSU-RAY1NC/SUTOM  1-25.6 

C 
C 
C 

C  READ  IN  MAP  ARRAY  DATA  FROM  DISK  FILE 
C 

C  BLKL4  DISC  ADDRESS  OF  HIPSOGRAPHIC  DATA 

0029  BLKL4-*056730 

C  MPBYTS  NO  OF  BYTES  OF  HIPSO  DATA  ON  THE  DISC 

0030  MPBVTS-32768 

0031  CALL  UDMAP(MAP,MPBYTS,BLKL4,IER) 

C 

C 

C 

C  FIND  MAXIMUM  HEIGHT  OF  TERRAIN  IN  MAP  ARRAY 

C 

C 

0032  HMAX  -  0 

0033  DO  10  J  -  l.NMAP 

0034  DO  10  1  -  l.NMAP 

0035  IF  (MAP(O.I)  .GT.  HMAX)  HMAX  ■  MAP(O.I) 

0036  10  CONTINUE 

0037  TYPE  MAX  HEIGHT  IN  MAP  ARRAY  IS: HMAX, MPBYTS, BLKL4, IER 

C 
C 

C  /***••  OUTER  PROGRAM  LOOP  **•**/ 

C  DO  ONCE  FOR  EACH  16  BY  16  KM  AREA 

C 

C 

0038  169  TYPE  INPUT  A/C  X.Y.ALT  (IN  UNITS  OF  KM  FROM  ALPHA, BETA)' 

0039  ACCEPT  *,XACM, YACM.ZACM 

0040  XAC-XACMM000.0/SUTOM  (CONVERT  TO  SCOPE  UNITS 

0041  YAC- YACM* 1 000 . 0/SUTOM 

0042  ZAC-ZACM* 1000.0  ICONVERT  ALTITUDE  TO  METERS 

C 

C  INITIALIZE  STARTING  COORDINATES  OF  SCAN  TO  LOWER  LEFT  CORNER 


13:27:00 
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FORTRAN  IV-PLUS  V02-51 
T.FTN  /WR 

0043  IX-IX1-IDEL 

C 

0044  INVISs0 

C 

C  THIS  IS  THE  MAIN  PROGRAM  LOOP  WHICH  ITERATES  ON  THE  NUMBER 

C  OF  VERTICAL  SCANS  (ROADS)  THAT  ARE  TO  BE  SEARCHED 

C 
C 
C 

004S  300  DO  8  NRD-l.NINV 

C  GET  STARTING  COORDINATES  FOR  BOTTOM  OF  NEXT  SCAN 
0046  IX-IX+IDEL 

0047  IY-IY1 

0040  DO  87  JK-1,65  iCLEAR  COLUMNS  1  TO  130  IN  MARK  ARRAY 

0049  87  1MARK< JK>»*20040  (WITH  ASCII  SPACES 

C 
C 

C  THIS  IS  THE  MIDDLE  LOOP  WHICH  ITERATES  ON  THE  NUMBER  OF  POINTS 

C  ALONG  THE  GIVEN  ROAD  THAT  ARE  TO  BE  SIGHTED  BY  THE  AIRCRAFT 

C 
C 

0050  DO  7  ITSDUM-3.NINV+2 

C 
C 

C  CONVERT  IX, IY  INTO  UNITS  OF  MAP  MATRIX,  AND  GET  HEIGHT  OF  ROAD  AT 

C  THIS  POINT  IN  CONTOUR  UNITS 

C 

0051  301  Xl-< IX+1024 )*XMAPS*1. 5 

0052  IXM-X1 

0053  Yl°( IY+1024 )*XMAPS+ 1 . 5 

0054  I YMBY 1 

CHECK  TO  SEE  IF  STARTING  POINT  IS  IN  MAP  ARRAY 
C  ASSUME  VISIBILITY  IF  NOT 

0055  IF( IXM.LE.0.OR.IXM.GT. NMAP )  GO  TO  6 

0056  IF< IYM.LE.0.OR. IVM.GT.NMAP)  GO  TO  6 

C 

0057  302  IIO«MAP<  IXM.IYM) 

0058  RA  =  RE  ♦  HO*CONTUR 

0059  HO  -  HO  -  TREE 

C 

C  COMPUTE  DISTANCE  FROM  A/C  TO  POINT  ON  GROUND  IN  SCOPE  UNITS 
C 

0060  303  DX-XAC-IX 

0061  DV-YAC-IY 

0062  DIST»SQRT(OX*DX+DY«DY) 

0063  IF(DIST.LE.0.0)  GO  TO  6  IAVOID  DIVIDE  BY  0  WHEN  A/C  IS  OVER  THE  POINT 

C 

C  GET  INCREMENTS  <  D!( ,  DY )  IN  MAP  UNITS  TO  FOLLOW  RAY  FROM  GROUND  TO  A/C 
0064  P0RT=XMAPS*R1NCSU/DIST 

0065  DX«DX*PORT 

0066  DY-DY*PORT 

C 

C  CONVERT  DISTANCE  TO  METERS  FROM  SCOPE  UNITS 
C 

0067  DIST  -  DIST  *  SUTOM 

C 

C  RA2  ANO  A  ARE  FUDGE  FACTORS  USED  TO  TAKE  INTO  ACCOUNT 
C  THE  FACT  THAT  THE  EARTH  IS  SPHERICAL 

C  RA2  IS  IN  METERS* *2/C0NTUR  UNIT 

C  A  IS  IN  METERS 
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FORTRAN  1V-FLUS  V02-81 
T.FTN  /WR 


13:27 >*R  02-OCT-79 


FACE  4 


C 

0069  RA2  •  2  *  RA  •  CONTUR 

#069  A  -  -  MOIST  •  OIST )  ♦  IRA  •  RA)  -  (RE  *  ZAC)  2)  /  OIST 

C 

C  INITIALIZE  VARIABLES  TO  BE  USED  IN  INNER  LOOP 
#07#  YINC  «  RAY INC 

C 
C 

C  INNERMOST  LOOP  STARTS  HERE 

C  THIS  LOOP  TRACES  A  LINE  FROM  OUR  POSITION  TO  THE  AIRCRAFT 

C  DETERMINING  IF  THERE  IS  ANY  TERRAIN  BLOCKING  THE  LINE-OF-SIGHT 

C 
C 

C  CET  COORDINATES  OF  NEW  POINT  ON  RAY  TO  AIRCRAFT  IN  MAP  UNITS 
C 

0071  400  X1“X1*DX 
0072  V1»YI*0Y 

0073  IXM  •  XI 

0074  IVM  -  VI 

C 

CHECK  TO  SEE  IF  RAY  IS  OUT  OF  BOUNDS 

007 S  IF( IXM. LE .0.OR . IXM.GT.NMAP )  CO  TO  6  IWE  HAVE  MOVED  ALONG  THE  RAY  SO  WE  ARE 

0076  IF( I VM.LE .0.OR. I VM.GT. NMAP  >  CO  TO  6  IOUTSIDE  THE  HIPSO  DATA.  CALL  VISIBLE 

C 

C  CET  HEIGHT  OF  WHERE  WE  ARE  ON  RAY  <Z1>  AND  COMPARE  WITH  HEIGHT 

C  OF  TERRAIN  AT  THIS  POINT 

C  21  IS  IN  CONTUR  UNITS  TO  SPEED  UPCALCULATION. 

C  A  AND  YINC  ARE  IN  METERS 

C 

0077  401  Zl“HO*Y JNC*( A*Y1NC )/RA2 

0078  VINC>YINC*RAV1NC 

0079  IF(Zl.GT.HMAX)  GO  TO  B  1D1D  ENOUGH.  THE  POINT  IS  VISIBLE 

0080  IF1Z1.GT. MAPI IXM, IYM)>  GO  TO  400  IMOVE  TO  NEXT  POINT  UP  THE  RAY  ANO  REPEAT 

C 

C  IF  WE  FALL  THROUGH  TO  HERE  THE  ROAD  SITE  IS  INDEED  BLOCKED 
C  FROM  THE  VIEW  OF  THE  AIRCRAFT  I 

C 
C 

0081  102  INVIS  -  INViS  ♦  1 

0082  MARKUTSDUM>-42  I  ASTERISK 

0083  IY*IY*IDEL 

0084  GO  TO  7 

C 

C  BUMP  ROAD  COORDINATES  BY  AN  INTERVAL  AND  LOOP 
C 

0086  6  IY-IV-HDEL 

0086  MARK< ITSOUMI-46  I  PERIOO  IF  VISIBLE 

0087  7  CONTINUE 

C 

c 

c 

0088  101  1F(LUN.NE.0)  WRITE  (LUN.234)  MARK  (OUTPUT  THE  RESULTS  OF  THIS  VERTICAL  SCAN 

0089  234  FORMAT! IX, 1 30A1 > 

C 

C  END  OF  ROAD  LOOP.  GET  NEXT  ROAD 
C 

0090  8  CONTINUE 

C 

C  END  OF  LOOP  OUTPUT  THE  RESULTS 

NV1S>100.*INVIS/FLOAT<MINV*N1NV)+0.S 
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0302 

0003 

0034 

000  5 


0096 


0097 


C 

1F(L«N.EQ.0>  TYPE  236,NVIS,MACM,VACM,2ACM 
236  FORMAT  (15,'  PERCENT  IS  INVISIBLE  ATs‘3F10.1> 

I F  ( L UN .  NE  .  B  >  V/R  l TE  <  MIN ,  233  )  XACi  1 ,  ‘'ACM .  '.’AIM ,  NVI 3 
235  FORMAT! '  AIRCRAFT  AT  •  .2F5.1,*  KM  FROM  CEMTER , ' ,F5. I , 
1  1  KM  ALTITUDE.1 ,14,*  X  IS  IN71SI3LE  *  > 

C 

GO  TO  169 
C 
C 
C 

END 
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DEC 
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GPC 

MASR 

POP 
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Digital  Graphics  Recorder 
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Queued  Input  Output 
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